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APPLICATIONS TO DETERMINISTIC AND STOCHASTIC 
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Abstract. In this paper we use a splitting technique to develop new multiscale basis functions 
for the multiscale finite element method (MsFEM). The multiscale basis functions are iteratively 
generated using a Green's kernel. The Green's kernel is based on the first differential operator of 
the splitting. The proposed MsFEM is applied to deterministic elliptic equations and stochastic 
elliptic equations, and we show that the proposed MsFEM can considerably reduce the dimension 
j^JQf of the random parameter space for stochastic problems. By combining the method with sparse grid 

i"^ ' collocation methods, the need for a prohibitive number of deterministic solves is alleviated. We rigor- 

ously analyze the convergence of the proposed method for both deterministic and stochastic elliptic 
equations. Computational complexity discussions are also offered to supplement the convergence 
analysis. A number of numerical results are presented to confirm the theoretical findings. 
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1. Introduction. Many fundamental and practical scientific problems involve 
a wide range of length scales. Typical examples may include subsurface flows and 
geophysical domains with microscopic structures. Because there exist both natural 
randomness and lack of knowledge about the physical properties, it is often necessary 
to incorporate uncertainties into the model inputs. One way to address the uncer- 
tainties is to model the random inputs as a random field/process, and in turn, such 
problems are often modeled as stochastic partial differential equations (SPDEs). Then 
the model's output can be accurately predicted by efficiently solving the associated 
SPDEs. It is challenging to solve the SPDEs when the random inputs vary over mul- 
tiple scales in space and contain inherent uncertainties. The interest in developing 
stochastic multiscale methods for the SPDEs has steadily grown in recent years (see 
e.g., 0[T3[ll[2TJ[26]). 

Let fi be a set of outcomes and D be a bounded domain in R d with a Lipschitz 
boundary. We consider the stochastic elliptic boundary value problem: seek a random 
field u(x, uj) : D x SI — > M such that u(x, uj) almost surely (a.s) satisfies the following 
equation 

( -V ■ (k(x,oj)Vu(x,oj)) = f(x) in D 
u{x,ui) =0 on 3D, 

where k{x,uj) is a scalar random field. In particular, we assume that k(x,oj) exhibits 
heterogeneity in multiple scales over space. Since k{x, uj) varies over different spa- 
tial scales, resolving the finest scale is not computationally feasible. Thus, we use 
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multiscalc methods. In practice, a high dimensional random field can be used to ap- 
proximate the stochastic field k(x, ui), yet computing the statistical output quantities 
of interest remains a difficult task. 

During the last decade several multiscalc methods have been developed, see e.g. [TJ 
[21 |3j OH [TUJ [15l [TT1 [18] . The idea of multiscale methods is to divide the fine scale field 
into many local sub-problems and solve these in order to form a global coarse scale 
equation. This leads to a coarse scale equation in which the fine scale effects are 
taken into account. One such multiscale method is the Multiscale Finite Element 
Method (MsFEM) [T5]. The main idea of MsFEM is to incorporate the small-scale 
information into the finite element basis functions and capture their effects on the large 
scale through the discrete variational formulation. In many cases, the multiscale basis 
functions can be pre-computed and used repeatedly in subsequent computations with 
different source terms, boundary conditions and even modified coefficients. 

The goal in this paper is to quantify the uncertainty through computing the sta- 
tistical moments (e.g., expectation and variance) of the stochastic solution. We note 
that the variance of the solution gives a measure for confidence of the solution expec- 
tation. Numerical solution strategies for stochastic PDEs generally follow three main 
steps. First, the random inputs are approximately parameterized by a finite number 
of random variables. This can be achieved by a truncated Karhuncn-Loevc expansion 
and/or truncated polynomial chaos expansion [25) . Second, a numerical approxima- 
tion for the resulting high-dimensional deterministic PDE is used to approximate the 
solution with the respective input parameters. Finally, the solution is reconstructed 
as a random field and the statistical quantities of interest are computed. The second 
step is most difficult because the PDEs involve high-dimensional random parameter 
inputs. There exist many methods for the second step. A broad survey of these 
methods can be found in [HI US]. Among these methods, Monte Carlo and stochastic 
collocation methods have been extensively studied and widely used. Monte Carlo 
methods and stochastic collocation methods generate completely decoupled systems, 
each of which is the same size as the deterministic system. This is suitable for parallel 
computing and amendable for relatively high-dimensional random inputs. In a Monte 
Carlo context, a large number of samples are randomly chosen and separate solves 
for each of the samples are used to determine the statistical behavior of solutions. 
However, a limitation is that convergence of Monte Carlo methods is usually slow. 
Unlike Monte Carlo methods, stochastic collocation requires independent solves at 
fixed collocation points which are specifically chosen. In turn, this type of method 
has the capability to provide better accuracy than Monte Carlo with a fewer number 
of realizations. Moreover, to overcome the curse of dimensionality imposed by high- 
dimensional input parameters, we can use Smolyak sparse grids (see e.g., 01231 [24]) to 
reduce the number of collocation points. In this paper, we consider the Monte Carlo 
method and the Smolyak sparse grid collocation method for stochastic approximation. 

In the paper, we consider both multiscale features and uncertainties simultane- 
ously A main focus is the use of a resourceful splitting technique to compute MsFEM 
basis functions. For the problem we assume that the coefficient k can be split 

into two parts, k = ko + ki. We then construct Green's functions using the differential 
operator associated with ko- The Green's functions are used to construct a sequence 
of multiscale "bubble functions," which are employed to build the multiscale basis 
functions for MsFEM in an iterative manner. The Green's function technique pro- 
vides an modified framework to compute the bubble functions, and is suitable for 
parallel computing due to the independent construction. The splitting of k is flex- 
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ible and can be easily controlled to lead to fast convergence of the bubble function 
sequence. Compared to standard MsFEM [15], the proposed MsFEM approach can 
accurately approximate multiscalc solutions. The new multiscale approach is applied 
to SPDEs and may result in new stochastic multiscale methods. Since the Green's 
functions essentially generate the MsFEM basis functions, this will reduce the dimen- 
sion of random parameter space if the dimension of the random field for ko is smaller 
than that of the random field for k. We note that using Karhunen-Loeve expansion 
or polynomial chaos expansion usually yields an inherent splitting. The new MsFEM 
can efficiently solve SPDEs with high-dimensional input parameters, and combining 
the approach with sparse grid collocation methods alleviates the need for a prohibitive 
number of deterministic solves. We present convergence analysis of the proposed Ms- 
FEM approach for deterministic elliptic equations and stochastic elliptic equations. 
Complexity analysis is also presented for deterministic MsFEM basis functions and 
stochastic MsFEM basis functions. 

The rest of the paper is organized as follows. In Section 2 we present the splitting 
technique which is used to compute the new MsFEM basis functions for deterministic 
elliptic PDEs, and provide the associated computational algorithm. In Section 3, 
convergence analysis is rigorously derived for deterministic elliptic PDEs. Section 4 
is devoted to the applications to stochastic elliptic PDEs. We present convergence 
analysis and complexity analysis using stochastic collocation methods in the section. 
In Section 5, a number of numerical examples are presented to confirm the theoretical 
results. Some conclusions and closing remarks are made in Section 6. 

2. A new approach for MsFEM basis functions. We consider the deter- 
ministic elliptic equation 

-V ■ (fcVii) = / in D 

V ' (2.1) 
u = on dD, 

where k is a heterogenous scalar function. We note that our method can immediately 
be extended to the case of tensor coefficient function. We assume that k admits the 
splitting, 

k = k + ki, (2.2) 

where k{x) and k$(x) arc bounded below and above, specifically, 

< a < k(x) < ai, < b < k a (x) < bi Vie D. 

Here, ko often represents the coarse scale information of fc, and k\ the fine scale 
information of k. 

The multiscale finite element method (MsFEM) for Eq. (|2.ip was introduced in 
[15] and further analyzed in [TB] . The key ingredient of MsFEM is the construction of 
an appropriate multiscale finite dimensional space in which the solution is sought. In 
particular, the fine scale heterogeneity in k will be imbedded in this finite dimensional 
space. This information is incorporated into the coarse scale formulation through the 
coarse scale stiffness matrix. In this section, we develop a MsFEM basis function, 
which is constructed in a different way from the previous works (e.g., [4] 115]). 

We introduce some notation for presentation. L P (D) (1 < p < oo) denotes the 
Lebesgue space. The norm of L 2 (D) is denoted by || • ||o,d- ^(D) is the usual 
Sobolev space equipped with norm || ■ \\x,d and scminorm | ■ \i,d- In the paper, (■, ■) 
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is the usual L 2 inner product. We define an energy norm on a sub-domain D' by 
IIMIId' : ~ (kVv,Vv) D > = ||\/feVi;||o D i- If D' = D, then ||| • ||| simply represents 
1 1 1 • 1 1 Id ■ We let Hh be a quasi-uniform partition of £l and K be a representative coarse 
mesh with diam(iT) = Kk- Let h — maxj/i^, K G 2^}. 

2.1. Series approximation of multiscale basis functions. Following [T3] , 
we define the standard multiscale basis functions by for vertices i = 1,. . .,d of 
coarse cell if, which satisfy 

r-V-(W^)=0 in if 

| <?f>KT,i = on dK, 

where Ik,i is the boundary condition associated with node i. There exist some options 
for the boundary condition Ik,% (see [HI [Til [20] ) . Eq. (|2.3|) defines basis functions for 
local MsFEM if Ika is a linear/bilinear function. Incorporating global information 
into produces global MsFEM [20] . We define the finite element space for the 
standard MsFEM by 

V h = spa,n{(/) K j : i = 1, d; K G 1 h }- 

We note that the idea of using basis functions satisfying certain differential equations 
has been used before, see e.g. [4] [19] and the references therein. Since we discuss 
a generic multiscale basis function, hereafter, we will remove the subindex K and i 
from (|2.3p for simplicity of presentation. 

Next we use the splitting (|2.2|) to derive a new MsFEM basis function. On each 
coarse cell K <G 2/,., we define a projection operator LI : iJ 1 (i4T) — > Hq(K) by 

(fcoVILv, Vw) = (fc Vw, Vw) Vv G H X {K) and Vw G H^(K). (2.4) 

The definition of IL implies |j\/fcoVLb;||o,if < || v^oVullo.if ■ 

We extend the function for the boundary condition in (|2.3[) onto K and denote 
the extended function by /. Let </>=(! — LI)/ + £, where I is the identity operator. 
Then by (|2.3|) we can derive an equation for £ 

f -V • = -V • (fciV(II - 7)0 in K 

\ ^ = on 3Jf. 1 ' j 

We arc going to construct a series to approximate £. To this end, we set £o to satisfy 
f-V-(AoV6) = -V-(AiV(n-I)0 in ^ 

\ f = o on air. 

Then we recursively define a sequence of function £j, j = 1, 2, 3, ■ ■ ■ which satisfies 

f-V-(fc V| i ) = V-(fciV| j _ 1 ) in If 
\ = on dK. 

The function LI/ and the sequence are "bubble functions" containing microstruc- 
ture information, which are localized to a coarse cell by imposing zero Dirichlct bound- 
ary conditions. Let £j = X)/=o£r We d erme the new multiscale basis function 
(pj := (I — LI)/ + £j and the finite element space for the new MsFEM by 

Vj. h = spa,n{(<j)j) K .i :i = l, d; K G T^}. 
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2.2. Computational approach for the new MsFEM basis functions. 

Since the proposed multiscale basis function is denned as <j>j = (I — H)l + Ylj=o£j' 
the computation of 4>j depends on the construction of ITZ and the bubble sequence 

{lj}/=o- We nnd that m and £3 (J = °> • • ■ > J ) 

are all associated with the differential 
operator £o := —V • fcoV. Moreover, ITZ or £j (j = 0, • • • , J) can be formally written 
as Cq 1 /, where / is the source term in the equation on HI or £j. It is well-known 
that Green's function can be viewed as generalized inverses of differential operators. 
We use Green's functions to obtain Cq 1 in the present work. 

As far as making an efficient implementation, we use the Green's function G(x, y) 
for the operator Co. The Green's function G(x,y) solves the equation 

f -V ■ (koVG(x,y)) = S(x,y) in K 

\ G(x, y) = on dK. ( ' ' 

Since the Green's function G(x, y) offers the fundamental solution for the differential 
operator Co, the Green's function G(x,y) can efficiently generate HI and (j = 
0, ■ ■ • , J). 

By Eq. ([23]), we have 

IB(x) = - / G(a;,y)V y • (fc V y /)dy = / fc V a G(x,2/) • V„Wy. (2.9) 

Then we similarly compute £o by performing 
io(x) = - f G(x, y)V v ■ (hV y (U - I)l)dy = f hV y G(x, y) ■ (V„(II - I)l)dy, 

J K J K 

and compute k = 1, 2, 3, • • • , by performing 

ij( x ) = / G(x,y)V y ■ (kiVyij-i)dy = kiV y G(x,y) ■ V v £j-idy. 

JK JK 

To discuss the complexity of computation for the proposed MsFEM basis function, 
we investigate the computation in terms of matrix operations. We use the vector 
function b(x) = (£i(x), • • • , £ nK (x)) T , where £ p (x) (p — 1, • • • Uk) is a standard finite 
element basis function at the underlying vertex x v of a underlying hne grid in K. 
Here tik is the number of the internal hne vertices in K. We define vectors vq and v\ 

by 

v = k Vb <g> Vldx and v x = / fc x V6 <8> VWx, (2-10) 

where (X) represents the tensor product. We use C\ = —V • fciV to denote the differen- 
tial operator associated with k\ . Let Mo and Mi be the stiffness matrices associated 
with the operators Cq and C\, respectively. Then 



M = / fc V6 ® VMi and M x = / kxVb® X7bdx. (2.11) 

We have the following theorem to represent 11/ and £j (j = 0, ■ • • , J) in the finite 
element space of fine grid. The notations 11/ and £j arc slightly abused in the following 
theorem. 
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Theorem 2.1. Let Hl(x) and £j(x) (j — 0, • ■ • , J) be the finite element approxi- 
mations on the underlying fine grid in K . Then 

IU(x) = (M " 1 6(x)) T u (2.12) 

and for j = 0, • ■ ■ , J, 

ij(x) = (-l) 3 '(M - 1 6(o:)) T (MxM - 1 y(MiM - 1 «o - (2.13) 



Proof. Let us still use G(x,y) to represent the numerical Green's function on the 
underlying fine grid. Then direct calculation implies that 

G(x,y) = (M^b(x)) T b(y). (2.14) 
Thanks to Eq. (gH) and Eq. ([2~T4]) . it follows that 

II/(ar)= / k V y G(x,y)-V y ldy= [ k^ y {M ( ; 1 b{x)) T b{y) -V y ldy 

Jk Jk (2.15) 

= (Mo 1 b{x)) T / k V y b®V y ldy = (M - 1 6(x)) T ?;o. 
■/if 

For £oOe)j we have 

£o(:c) - / hV y G(x,y) ■ V y ni{y)dy- [ kiV y G(x,y) ■ Vyldy 
Jk Jk 

= [ k 1 Vy{M^ 1 b{x)) T b{y)-Vy(M^ 1 b{y)) T v dy-{M^ 1 b(x)) T j k x V ' y b(y) ®V ' y ld 
Jk Jk (2.16 

= (M^b{x)) T [ f hV y b(y) ® V„? r (y)dtf]M - 1 «o - (M^ 1 b(x)) T v 1 
Jk 

= (M - 1 6(x)) T (M 1 M - 1 v - 
Using Eq. (|2.16p . we obtain 

£ 1 (x) = - kiV y G(x,y) -V y Co(y)dy 
Jk 

= -(A/ - 1 6(x)) T [ / hVybiy) ® V y b r (y)dy]Mo 1 (M 1 Mo 1 v - v^ 211 ) 
Jk 

= -(M^b{x)) T {M 1 Mo l ){M 1 M^ 1 v - Vl ). 
By repeating the procedure of (|2.17[) , it follows immediately that for j = 2, • • • , J, 

lj(x) = (-l) 3 '(M - 1 6(o;)) T (M 1 M - 1 y(M 1 M - 1 i;o - ui). 

The proof is complete. □ 

We pre-compute vectors Vq, i>i and matrices Mo, Mi. Since vq, v\, Mq and M\ 
only depend on the local information in K, their construction is suitable for parallel 
computation. By Theorem l2.ll we obtain Tll(x) and £j(x) (j = 0, • • • , J) by perform- 
ing a direct matrix- vector multiplication. Moreover, we find that the computations of 
Hl(x) and £j(x) (j = 0, ■ • • , J) are independent of each other and suitable for parallel 
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computation as well. By Theorem 12. 1[ the numerical representation of <j)j can be 
written as 

j 

<f>j(x) = l(x) - (M - 1 ?(x)) T «d + J2(- 1 ) j ( M o 1 b(x)) T (M 1 Mo 1 ) j (M 1 Mo 1 v - Vl ). 

3=0 

(2.18) 

If we still use 0(x) to denote the numerical approximation of a standard MsPEM basis 
function in the underlying fine grid of K, then it is easy to show that 

4>(x) = l(x) - [M~ l b{x)) T v, (2.19) 

where M = J K kVb® Vbdx and v = J K kVb® Vldx. Compared Eq. (|2.18p and Eq. 
(|2.19j) , we find that the computation of 4>j{x) is comparable to the computation of </>(x) 
in a parallel setting since each of the terms in 4>j(x) may be obtained independently. 
When ko = ki, Eq. (|2.18|) and Eq. (|2.19p imply that <j>j{x) = 4>{x), which means 
that the proposed MsFEM coincides with standard MsFEM. 
Application of this procedure is summarized in Algorithm [T] 



Algorithm 1 Computing the proposed MsFEM basis function <fij 

1. Assemble the vectors vo and v\ by (|2.10[) . and assemble the matrices Mq and 
Mi by Q3D; 

2. Compute each term of ()2.18j) independently; 

3. Construct the proposed MsFEM basis function 4>j(x) by summing the terms 
obtained in step (2). 



Remark 2.1. Let £ be the finite element solution (i.e., evaluated at fine vertices) 
of Eq. h2.5\) . Here the notation £ is slightly abused. Then by Eq. \2. 5)) we have 
(Mq + Mi)£ = MiMq 1 vq — Vi. This gives the Neumann series of 

oo 

t = E^W^iMT WW^o - vi). (2.20) 

3=0 

Consequently, the last term in Eq. i2.18\) is actually the truncated Neumann series. 
By the Neumann series and Eq. i2.18\) . if the spectral radius of (MiMq 1 ) is less than 
1, then <pj is convergent as J — > oo. This will be consistent with Assumption 3.1. 

3. Convergence analysis for deterministic coefficient . We make the fol- 
lowing assumption for convergence analysis. 

Assumption 3.1. Assume that the splitting k = ko + ki on K satisfies 

r)K-=\\^\\L-{K) < 1. (3.1) 

If the property of the splitting ||!^||l°°(j<-) < 1 fails, then we can introduce a slightly 
modified splitting of k such that the assumption (|3 . 1 1) is satisfied. As a matter of fact, 
we have the following proposition. 

PROPOSITION 3.1. Suppose that \\jr\\L<*>(K) > 1- Let constant s satisfy the 
following inequality 

( k\ ix) — kr\(x) , , A . . 

s > sup max n , ~k (x) . (3.2) 
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Then the modified splitting k — fc + k\ := (k + s) + (fcj — s) satisfies the assumption 
1)3.1)) if we switch (fco,fci) to (ko,ki) in \3.1)) . 

Proof. If H^Hioo^) > 1 and (|3.2p holds, then direct calculation implies that 

II fHU^tA') < 1- Consequently, the assumption (|3.1[) is satisfied. □ 

We note that the modified splitting, k = k + fc\ , is essentially the same (up to a 
constant) as the original splitting, k = fco + fci . 

Now we provide the main convergence result for deterministic MsFEMs. 

Theorem 3.2. Define c = \f2max-K&t h \\\J^;\\l'={k) andn = maxxei h Vk < 1- 
Let Uh G Vh and uj t h G Vj,h be the MsFEM solution for Eq. 112.1)) . Then 

(1) UK - Uj, h \\\ <(o)^ 2 {v J+1 + T ^) 1/2 ||| U |||. 

(2) \\\u-uj, h \\\ < V2£(r] J+1 + 1 * J v + , 1 +1 )\\\u\\\+2\\\u-u h \\\. 

Thanks to the boundedness of fco and fc, it follows that c < By Theorem 

it follows that 



lim \\\uh — Ujh\\\ = 0, lim lim |||u — = 0. 

J— loo h— >0 J— too 

This shows convergence of the proposed MsFEM. 

The arguments of the proof Theorem [3~2l include the following Lemma l3.6l Lemma 
3.71 and Lemma 13.81 In the rest of the section, we will formulate these lemmas. 

To describe and prove Lemma 13.61 we need a couple of lemmas. The following 
lemma gives an upper bound of the sequence {£■/}• 

Lemma 3.3. Let £o and £j (j = 1, 2, • • • ) solve Eqs. \2. 6)) and {2.7)) , respectively. 
Then 

feV| J -|| ,K<2< +1 ||V^u'VZ||o,A, j = 0,1,2,-.. , 



where tjk is defined in \3.1)) . 

The proof of Lemma 13.31 is presented in Appendix El 

The following lemma shows that the sequence of the new MsFE basis functions 
converge to a standard MsFE basis function. 

Lemma 3.4. The sequence of basis functions {4>j} is convergent. Moreover, 

lim 111^-011^ = 0, 

J— foo 



where </> solves the standard MsFE basis equation S2.3)) . 
The proof of Lemma 13.41 is presented in Appendix [EJ 

From the proof of Lemma l3.4| we can obtain an explicit convergence rate for <pj, 
which is stated as following: 

Proposition 3.5. If the coarse function I in Eq. \2. 6)) is in H 1 ^), then 

|||0-^||k<C^^+ 2 , 



where the constant C'i = ||V/||o,a". 

The proof of Proposition 13.51 is presented in Appendix O 

To describe Lemma [3~51 we need to introduce elliptic projection operators. Specif- 
ically, we define the elliptic projection Hh : Hq(D) — > Vh by 

(fcVn h v, Vw) := {kVv, Vw) W G H%(D) and Vtu G V h . 
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Similarly, the elliptic projection Hjh '■ Hq(D) — > Vjh is defined by 

(kVIlj >h v, S7w) := (kVv, Vto) Vu G iTo(-D) and Vw £ V J)h . 

Since limj^oo Vj t h = Vh by Lemma 13.41 the standard density argument [8] implies 
that limj-j-oo Tlj.h — Tih- We can show that the projection operators Tlh and Tlj,h are 
self-adjoint and idempotent operators. Let uk and uj,h be the MsFEM solution in Vh 
and Vj,h, respectively, for Eq. (|2.1|) . Then 



u h = U h u and uj, h = ILj ih u. 

Using Lemma [3.3l and the technique presented in [14] . we have the following result. 
Lemma 3.6. Letue H^(D). Then 

iij^iMII + Ilia - n h )n Jlfc u||| < c(n J+1 + ^_ + j +1 )|||u|||, (3.3) 

where c and rj are defined in Theorem \3.2i 
Proof. For any u G Hq(D), we write 

K i K i 

We define the operator T : H%(K) — > H%(K) by 

-V • {k VTv) = V • (fciVw), y-v G #o(-*Q- 

Then we have 

oo oo 

<l>K,i-{<i>j)K,i= 51 (C?W = TJ (£o)if,i 
j=J+l J=J+1 



(3.4) 



T j+1 J2 T3 (£o)k^ = T J+1 {4> K ,i -{I- T\)l K .i). 

3=0 



By (|3.4p . we have 

(7 - nj^lMx = ^ a K ,i {<j> Kti - {<f>j) K ,i) = T J+1 (n h u\ K - £ OkA 1 ~ U ) l K,i) ■ 

i i 

(3.5) 

and 

{I-Il h )Ilj, h u\ K = Y J PKA(<t>j)K,i-<t>K,i) = T ,+1 {Y.PKAl-n)l K , l ~U h Iljj l u\K)- 

i i 

(3.6) 

We first estimate (I — ilj^)n^u. Let Zk = ^Iiu\k — J2i a K,i(I — n)Z/r ». Thanks to 
T5|) and the proof of Lemma 13. 3[ it follows that 

(/ - n Jifc )n fc u|||^ = (kv(n h u - u Jth n hU ), v(n h u - n. Lh ii h u)) K 

= (kV(T J+1 z K ),V(Ii h u-Uj. h Ii h u)) K 

< ||V*v(T- /+1 2( J c)||o,jf|||(J-nj, fc )n h u||| J c 
T 



< \U Tr\\L^(K)\Wk Q V{T J + 1 Z K )\WK\W - ^j,h)^ h u\\\ K 



< v/2^ +1 ||Vfeu'v^||o 1 if|||(/-nj, /l )n /l u||| i f. 
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This implies that 

|||(7 -U Jih )IL h u\\\ K < V2 V i+ 1 \\^k~ Q Vz K \\ Q . K . (3.7) 
Since Zk € Hq(K) and (k S7(I — H)ljc t i,ZK) = 0, we have 

II VkoVz K \\l, K = (k U h u\ K - k ^2a Kii (I - H)Ik,uVzk) 

i 

= (k U h u\ K ,\/z K ) < ||y ^Ilioo^llinfculllxll-v/^oVarjfllo.jif, 



from which we get 



Hv^Vzjc||o,jc < lly ylli-(je)lll n h u llk- (3-8) 



By (|3.7I) and (|3.8[) . it follows immediately that 



(7-n Jlh )nAu||| < \/2maxM| A / 1 ^|| i =o (x) r/^ +1 J |||n fc tt| 
< CT7 J+1 |||u| 



(3.9) 



Next we estimate (I — HhjHj^u- By (|3.6|) . we have 

fe^v(/ - n fe )n J)hU || 0iK = ||v / ^VT J+1 (5^ J 8 J f ii (/-n)j jrii -n h n J)/i «| jr )||o 1 j C 

r (3.1-p) 

< Hv^v(nj, h u - 53#f,i(i - n)/ J f, i )||o,jc + ||v^v(n Jlh « - n^n^Uo,* 

Let z A ' = nj )fe M|if - £\ fe,j(7 - H)lic,i- Then f # e H^(K) and 

(fc V5if, Vzk) = (fcoVnj^w, Vzjf) < || V^oVIIjjjujlo.A'll V^oVzk||o,k- (3-11) 
Combining (|3.10[) and (|3.11|) implies that 

>— n J+1 j— 

i - 

Because (k VIlj !h u, VIL Jih u) < ||^||i«(jf)ll|n ./^ w lllifi ^ follows 

||V^v(J-n fc )n Jlfc «|| ,K < /Jf j +1 ||^IUco(joll|nj, h u|||A:- (3.12) 

1 — r\ K K 

Since |||^||z,~(K') < lj wc get || ^ |U°°(jr) > |- A direct calculation gives 

^v(i-n h )n Jifc ii||g )Jf > 5lH(/-n h )iij, h «|||^. (3.13) 



By (|3~T2|) and (|3~T3| . we have 



J+i 



(/-n fc )n Jih «|lk < v2||-rll^°°w -, j+i lW^hu 

1 - Vk 



K ■ 
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Since |||nj^w|||if < it follows immediately that 

IIKJ-n^nj^Hi < g 1 ^j +1 ||lu||k. (3.W) 

Combining p.9[) and p,14[) completes the proof. □ 

A straightforward application of Theorem 3.6 in |14j gives rise to the following 
lemma. 

Lemma 3.7. Let the inequality i3.3\) in Lemma \3.6\ hold. Then 



where c is defined in Lemma \3.b\ 

Lemma 13.71 gives the first convergence result in Theorem [ 

Following the proof of Theorem 3.8 in [T3] and Lemma |3"771 we have the following 
lemma. 

Lemma 3.8. Let u be the solution to Eq. h2.1}) and u h <G Vj, the standard MsFEM 
solution of Eq. h2.1\) . Then 

n J+1 

\\\u-TLj h u\\\< V2cU T+1 + — — )||| u ||| + 2|||u- u h \\\. 

1 — n J+L 

Because uj.h = Hj ,hU, the second convergence result in Theorem 13.21 follows 
Lemma 13.81 immediately. 

4. Analysis for stochastic coefficient. We consider the stochastic elliptic 
equation (jl.ip . To make k(x,ui) positive, we consider a logarithmic stochastic field, 
k(x,uj) := exp(Y(x, w)), where Y(x,ui) is a stochastic field with second moment. 

We assume that Y(x, uj) admits the following truncated Karhunen-Loeve expan- 
sion (see [25] for details), i.e., 



Y(x,u;)=E[Y] + ^2^ i bi(x)e i (ui). (4.1) 



i=i 

Here (h,bj) L 2( D -) = Ai > A 2 • • • , > A m • • • , lim m — >ao A m = 0. Let 9 := 

(01,-.. ,e m ,e m +i,--- ,0n) := (9o,ei) e R n , where 6 := (6i,--- ,9 m ) € M m and 
0i e M. n ~ m . Then stochastic field k(x, lu) can be parameterized to a finite-dimensional 
random field k(x, Q). We define the splitting of k(x, O) by 

*(x J 0) = Ao(x j e o ) + ft 1 (a:,e), (4.2) 

where k (x, @o) '■= ex P(-^[Y]+X)i!li y/\bi{x)0i{ijj)) (m < n) and fci(x, 6) = k(x, 6) — 
k (x, O ). 

4.1. Convergence analysis. In the subsection, we will present convergence 
analysis when the random field k(x, 0) admits the splitting described in (|4.2[) . 

We define L 2 (f2) to be the square integrable space with the probability mea- 
sure p(Q)dio, where p(0) is the joint probability density function of 0. We can use 
Theorem 13.21 to derive an error estimate of the stochastic elliptic equation (jl.ll) . 
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Theorem 4.1. Let and Ujh be the MsFEM solution of stochastic elliptic 
equation in space Vh x L 2 (il) and Vj.h x L 2 (fl), respectively, for the stochastic equation 

& Vv ■■= ll^gfjlU»(^xn) < 1, then 



\\VkV(uh - wj,h)||i2(£>)®L2(a) < {c) 1/2 {v J+1 + 1 _ J+1 ) ||VfcVu|| L 2 (£) ) (8L 2 (n ), 

^ ere c=V2||^||i|| L c. (Dx „). 

The eigenvalues {A^} play an important role to control To this end, we define 
the energy ratio E(m) by E(m) = 2 ~^i= 1 J^ . Then we can show that II ko(x,'&o) IU°°(°xn) 
is proportional to 1 — E(m) under certain conditions. 

PROPOSITION 4.2. Let ||^||l~(o) < C? uniformly for all m + 1 < i < n. If 
cov[Y](xi, X2) are piecewise analytic in D x D, then there exists constant a Cy such 
that for m large enough 

ll^b-^a-Ci-sM), 

where C Y = \Cg max m+ i< i <„{|6 i | i ^ (D )}. 

Proof. Since fci = fco(^ — 1) = fco(cxp(logfc — logfco) — l), then it follows that 

fc 



exp ^2 y/\ibi{x)9i(bS) - 1 



\i— m+1 



Because cov[F](a;i, x-i) is piecewise analytic in D x D, then there exist positive con- 
stants Co independent of m and <i (see [12]) such that for given < s < i 

n 

|| J2 V^Mx)Oi(u})\\Loo (Dx n)<C 1 exp(-C (--s)mi), (4.3) 



i— m+1 

where Ci = Ci(Cg,s, d, C ). If m is large enough, inequality (|4.3|) implies that 



^2 V^ibi(x)9i(uj))\\ L <x,( Dx a) < 1- 



i— m+1 

Due to the inequality le 3 ' — 1| < ||a;| for |x| < 1, then it follows 

||T-|U~(Dxn) < 7II \AA(z)0*H)|U°°(r>xn) 

i=m+l 

7 - 

< -C e max {Nl°° (£>)}( >Z V 7 ^" 

7 

< -C e max {|6 i | i oc (D) }(l-£;(m)). 

4 m+l<2<n 

The proof is completed. □ 



(4.4) 
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Theorem 14.11 and Proposition ^. 21 show that the convergence rate of the proposed 
MsFEM for the stochastic equation (jl.lj) depends on the energy ratio E(m). Then 
we immediately have the following theorem. 

Theorem 4.3. Suppose that the assumptions in Proposition \4-S\ hold. If m is 
large enough such that CV (l — E(m)) is less than 1, then 

\\VkV(u h - uj. h )\\ L 2 {D)(SL 2 {n) < C(C Y (1 - E{m)))^. 

If the stochastic field Y(x, oj) is Gaussian, then its covariance function can ana- 
lytically be extended to the whole complex plane C d , which is stronger than piecewise 
analytic. The eigenvalues {A^} associated with Gaussian fields decay very fast. Con- 
sequently, 1 — E{m) will decay very fast as m increases. Due to central limit theorem, 
the Gaussian stochastic process is very interesting for applications. 

For stochastic simulation, we can use Monte Carlo methods. The main disad- 
vantage of Monte Carlo methods is slow convergence. To overcome the disadvantage, 
here we use stochastic collocation methods to discrctizc random parameter space. 
The proposed MsFEM is used to discretize the spatial variable. Combined with the 
new MsFEM, we develop modified stochastic collocation methods to reduce the di- 
mension of the random parameter space. Computational complexity is addressed for 
the modified stochastic collocation method. 

4.2. Stochastic collocation methods and random parameter space re- 
duction. In this subsection, we use stochastic collocation methods to discretize the 
random parameter space and show that combining the proposed multiscale method 
with stochastic collocation methods can reduce the dimension of the random parame- 
ter space, which is very important for simulations in high-dimensional random space. 

Let {Oj, 8§, • • • , Qq} C R m be s collocation points scattered in random param- 
eter space associated with an interpolation operator I m . Let v(Qq) £ C(R m ) be a 
deterministic solution depending on random parameters Oq. Then given a realization 
0o £ R m , the collocation solution v m is defined by u m (©o) = I m v(®o). We usually 
use the roots of an orthogonal polynomial (e.g., Hcrmitc polynomial or Chebyshev 
polynomial) to find the interpolation points. One can select different collocation 
points and use a different interpolation operator I m to obtain different stochastic col- 
location methods, for example, full tensor product collocation [6] and Smolyak sparse 
grid collocation [TJ. 

Suppose that the Green's function G(x,y,Qo) m Q2.8P depends on the random 
parameter 0o £ R m and I m G(x,y,Qo) is the collocation solution for an arbitrary 
6 £ M. m . For any arbitrary realization 9 := (8 ,Qi) £ K m x W l ~ m , we define a 
modified interpolation operator I m for 11/ and £j (j = 0, • ■ ■ , J). Wc define them as 
follows: 

7 m [(ffl)(a;,eo)] := / k (y,G )V v I m G(x,y,G ) -W y ldy. (4.5) 

J K 

Wc note that wc can get the value of ko(y, 0o) and do not employ interpolation I m 
for ko(y, Qq)- Then we similarly compute 7 m [^o] by 

7 ro [lo(z,e)] := J h(y,Q)VyI m G(x )y> Q ) ■ (v v I m [QU)(y,Q )]-V v l(y)jdy, 
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and compute 7 TO [£fc], k = 1, 2, 3, • • ■ by performing 

J m [e/(a:,e)]:=- / fcifc,, e)V J/ 7 m G(a;,y, O ) • V^fe-ifo, 0)]d». (4.6) 

By the definitions of 7 m (H7] and I m [£,j} 0' = 0, • • • i 7), we have only used I m (G(x, y, 6o 
to compute 7 m [D7] and 7 m [£j] (j = 0, • • • , J). This interpolation is performed on 
the relatively low dimensional parameter space M m . Since the new multiscale basis 
function is defined as (f>j(x, 0) = i(a;) — Hl(x, 0g) + S/=o £:7'( a: > ®)> * nc computa- 
tion of the interpolation for (f>j(x, 0) only involves the m-dimensional interpolation 
I m (G(x,y,Q ))- 

To address the complexity, we use matrix and inner products to discuss the com- 
putation of the interpolations in (|4.5p and (|4.6[) . Due to (|2.12[) and (|4.5[) , it follows 
that 

7 m [m(x,0 o )] = A/^Mo-Heo)]^)) «b(e ). (4.7) 

By (j2~T3|) and gU), we have for j = 0, ■ ■ ■ , J, 

7 m fo(x, 9)] = f[/ m M - 1 (e )]6(s) > ) fM 1 (0)[7 m M o - 1 (0 o )]V 

(4.8) 

x f A/ 1 (e)[7 m Af- 1 (e )]wo(e ) -«i(e) 

By l|477) and (|4~8j) . we compute 7 m M o ~ 1 (0 o ), u (©o), «i(0) and Mi(0) to obtain 
7 m [(nZ)(x, ©o) and I m [£j(x, 9)] (j = 0, • • ■ , J). The computations of i>o(0 o ), i>i(0i) 
and Mi(0) are independent each other and very efficient in parallel. The dominant 
computation lies in 7 m M c j~ 1 (0o) and solely depends on the dimension m of the random 
parameter space. We can very efficiently compute 7 m M c j _1 (0o) in parallel as well. 

If we use the standard multiscale basis function defined in equation (|2.3[) . then 
the basis function is <f> :— (f>(x, 0) for an arbitrary realization £ R™. We interpolate 
the basis function (j>(x, 0) in the full random space IR ra (n > m). If n is large, the 
interpolation on M™ is computationally expensive and prohibitive. 

If we use full tensor product collocation and polynomials with degree q for each 
component of £ M n , then we have (q + l) n collocation points for full-space inter- 
polation. Consequently, we need to compute (q + 1)™ multiscale basis equations <j) for 
each vertex for the collocation. However, if we use the technique of the Green's func- 
tion for the new multiscale basis functions, then we have (q + l) m collocation points 
and compute only (q+ l) m Green's functions (or Green's matrix Mo) to generate the 
proposed multiscale basis functions. Consequently, the ratio a/t c between the number 
of collocation points in the proposed multiscale basis function and in the standard 
multiscale basis function is 

a ftc = ( q + l) m - n . 

For example, if n = 20, m = 10 and q = 2, then the ctftc = 3~ 10 . This means that the 
computational effort for using the proposed multiscale method may be considerably 
decreased compared to the standard multiscale method when the full tensor product 
collocation is employed. 
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Let H(n-\-L, n) denote the interpolation nodes for Smolyak sparse grid collocation 
[7] at dimension n and interpolation level L. Although Smolyak sparse grid collocation 
requires much fewer nodes than the full tensor product collocation to achieve the 
similar accuracy, the number of nodes H(n+L, n) increases very quickly as n increases. 
The ratio otftc between the number of collocation points in the proposed multiscale 
basis function and in the standard multiscale basis function is given by 

H(m + L,m) .m. T . 

a *9c = „( ; T — f w — for m > 1, n > 1, 
H(n + L,n) n 

where we have used the fact H(n + L,n) ps %jn L for n > 1 0, For example, if 
we take interpolation level L = 2, m — 10 and n = 20, then a sgc = ||j ps |. This 
means that the computation time of the proposed multiscale method is almost j of 
the standard multiscale method in parallel setting. 

Smolyak sparse grid collocation is known to have the same asymptotic accuracy as 
full tensor product collocation, while requiring many fewer interpolation points as the 
parameter dimension increases. We will use Smolyak sparse grid collocation for the 
numerical tests. The stochastic approximation of the Smolyak sparse grid collocation 
method, \\v — I m v\\ , depends on the total number of sparse grid collocation nodes and 
the dimension m of the random parameter space. The convergence analysis in |23j 
implies that the convergence of Smolyak sparse grid collocation with respect to the 
number of Smolyak nodes is exponential, but depends on the parameter dimension 
to. If to » 1, then the exponential convergence rate behaves algebraically. 

Using the modified stochastic collocation method described in (|4.5[) and (|4.6[) , we 
define the corresponding collocation basis function <j>j(x, 0) = l[x) — I m [Ul(x, ®o)] + 
^2j=o Im[ij(x, ©)]• Let uj.h be the collocation solution using the basis <j)j(x, 0). 
Then the total error 

— uj,h)\\L 2 (D)»L 2 (n) includes two parts: splitting error 

— u J.h)\\L 2 (D)®L 2 {n) and collocation error \\^fkV{u,j t h — ujji)\\l 2 {d)®l 2 (q,)- 
It can be formally expressed by 

— Uj,h)\\L 2 (D)(g>L 2 (n) < Cspl + e co i, 

where e sp i = 01 (l — E{mf) 2 I and e co ; = e co ;(J, m, L). There exists a trade-off 

between the splitting error e sp i and the collocation error e co ;. The numerical results 
in Section 5 illustrate the finding. 

Using the proposed stochastic collocation methods, computation of the MsFEM 
basis function <pj is summarized in Algorithm [5J 

5. Numerical Results. In this section we offer a number of representative nu- 
merical results to verify the analysis and evaluate the performance of the proposed 
method. 

5.1. Deterministic basis function. To begin, we consider the analysis offered 
in Sect. [3] In particular, we are initially interested in verifying the convergence prop- 
erties of a single, deterministic basis function as described in Lemma. 13.41 (or Prop. 
I3.5|) . In this subsection we consider two distinct cases of coefficient examples. We 
first consider a coefficient generated by a Karhunen-Loeve expansion and a coeffi- 
cient constructed from a log- normal distribution. See Fig. I5.1l for illustrations of each 
coefficient plotted on the log scale. 
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Algorithm 2 Computing MsFEM basis functions <j>j using the parameter reduction 
collocation method 

1. Choose collocation samples {Oo}f = i C IR m ; 

2. Assemble Green's matrices Mo(Oq) (i — 1, ■ ■ ■ ,s) for all collocation samples 
independently; 

3. Given an arbitrary realization := (©o, Oi) £ K m x K n_m , assemble vectors 
Wo(@o) and vi(©i) and matrix Mi(0) independently; 

4. Compute / m M ( ^" 1 (6o) indpendently; 

5. Compute I m [(W)(x,Q ) by (S3) and I m [^(x,Q)} (j = 0, ■ • • , J) by (jUJ 
independently; 

6. Construct the interpolated basis function <pj(x,Q) :— l(x) — I m [(Hl)(x, ©o) + 

E/= e)]. 




Fig. 5.1. ifLi? fie/tj and log-normal (right) coefficients posed on a 30 X 30 mesh (log scale) 



To generate the first test coefficient on an arbitrary coarse element K we employ 
the KLE expansion from Eq. (|4.ip . In our case we use the correlation function 

cov[y](x 1; x 2 ) := R{x x ,y x] x 2 ,y 2 ) = a 2 exp {- ^'^ - ^7^) . (5-1) 

where cr 2 is the variance, and l x ,ly denote the correlation lengths in the x— and 
y— directions, respectively We consider an elliptic coefficient which is generated on a 
30 x 30 mesh. For the variance we use a 2 = 2.25 and for the correlation lengths we 
use l x = 0.2 and l y = 0.05. For all examples we truncate the original KL expansion 
at n = 20 terms to obtain the full coefficient k = ko + ki . Then, in order to split the 
coefficient accordingly, we may choose a variety of m values and employ Eq. (|4.2j) for 
the splitting. For example, we may use m — 5 terms to obtain /co, and k\ — k — ko- 
See Fig. 15.21 for a representative example of a KLE coefficient splitting. As the initial 
analysis is built in a deterministic setting, we use the same, fixed 64 (i = 1, . . . , n) in 
(|4.2j) for all related examples. To begin, we recall the error estimate 

|||£-£/||k <n^\\L~(K)V J K +1 \\Vk~ Vl\\ , K (5.2) 
y/kkg 

and set := 2\\ ^ ||L° c (if) 7 7jf +1 ||v / fcoVZ||o,.K- from the last lines in the proof of 

Lemma. 13.41 In particular, we recall that when t\k = Wj^Wl^^k) < L convergence of 
the basis function sequence {4>j}y =0 is expected from the analysis. To test the error 
bound in Eq. (|5.2p we consider a variety of field splitting configurations. In Fig. 15.31 
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Fig. 5.2. KLE coefficient decomposition posed on a 30 X 30 mesh; n = 20, m = 5 




Fig. 5.3. Energy error and error bound computations for the KLE coefficient; m = 15 terms 
(left), m = 17 terms (right) 



we illustrate two cases of splitting where t]k = || I 1 ||l°o(js:) < 1- These examples result 
from the cases where m = 15 and m = 17 terms are used in the KLE splitting. First, it 
is important to note that the bounds presented in the analysis are clearly represented 
in the figure. In particular, for either case we see that the energy norm of the error 
is always bounded above by the theoretical estimate provided in Eq. (|5.2p . We also 
note that as J (the number of terms in the approximate basis function sequence) 
increases the error and associated bounds rapidly decrease. This behavior is expected 
as the term quickly decreases as J increases. We also point out that a smaller 

value of r/K yields a tighter bound. To conclude this coefficient example, we offer 
an illustration of the actual basis functions which are obtained through the proposed 
computational method. Fig. l5.4l includcs the benchmark basis function <fi as well as <pj 
(J = 0, 1, 2) from the sequence {(f>j}- We also plot the benchmark perturbation <f> — I 
and 4>j — I for J = 0, 1, 2. All plots were obtained from the case where m = 5 terms 
were used in the KLE splitting. We note that for a relatively pronounced splitting, 
we see a noticeable convergence trend. 

To generate a second test coefficient on an arbitrary coarse element K we assume 
that the full coefficient k follows a log-normal distribution. That is, we assume that 
]n(k(x)) = Y(x,w), where Y(x,uj) is a normal random variable with zero mean and 
variance one. We also assume that ln(fco(x)) = s c Y(x,u>), where s c is the strength 
factor which will determine the final splitting k = ko + fci. In particular, we set 



ki = k — ko = cxp(y) — exp(s c K) 



(5.3) 
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Fig. 5.5. Log-normal coefficient decomposition posed on a 30 X 30 mesh; s c = 0.4 



to create the coefficient decomposition. For the following examples, we choose various 
values of s c within the interval 0.4 < s c < 0.96. See Fig. 15.51 for an example of the 
coefficient splitting in Eq. (|5.3p for the case when s c = 0.4. We note that the fco 
portion of the decomposition is much less heterogeneous than k\. 

To validate the convergence properties outlined in Thm. 13.41 we consider two 
approaches and recall the estimate from Eq. (|5.2[) . In Fig. 15.61 we illustrate two 
representative plots of splitting where t\k = || j^\\l°°{K) < 1- These examples result 
from the cases where strength factors of s c = 0.90 and s c — 0.94 are used in the 
log-normal coefficient splitting. As before, we note that the bounds presented in the 
analysis arc clearly represented in the figure. For either case we see that the energy 
norm of the error is always bounded above by the theoretical estimate provided, and 
that as J increases, the error and associated bounds decrease rapidly. Here, we are 
also interested in the rate of convergence offered in Eq. (|5 . 2[) . To address this, we fix 
a value of J and plot t\k vs. \\\4> — 4>j\\\k on the log scale. Fig 15.71 illustrates the log 
plots as well as the slopes obtained from a linear trend line. In this case we obtain 
slopes which are close to the value of J + 2. In particular, for J = we obtain a 
slope of 1.8, for J = 2 we obtain a slope of 3.7 and for J = 4 we obtain a slope of 
5.6. These results are consistent with the convergence rate results from Prop. [3~51 In 
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Fig. 5.7. Log plot of r\ k vs. \\\<f> — <j>j\\\K for the log-normal coefficient; J = (left), J = 2 
(center), J = 4 (right) 



particular, we see from the plots in Fig. 15.71 that exponents of J + 2 — 5 are recovered 
for all values of J. Because the estimate in Prop. 13751 depends on a constant (maybe 
not sharp), there exist a slight difference S of the convergence rate in the numerical 
results compared to the estimate in in Prop. 13.51 

To conclude these coefficient examples, we offer a representative illustration of 
the actual basis functions which are obtained through the proposed computational 
method. Fig. 15.81 includes the benchmark basis function (j> as well as <f)j (J = 0,4,8) 
from the sequence {4>j}- We also plot the benchmark perturbation values 4> — I 
and 4>j — I for J = 0,4,8. All plots were obtained from the case where a strength 
factor of s c = 0.4 is used in the coefficient splitting. We note that this is a rather 
extreme splitting where fco is much less heterogenous compared to k\ (refer back to 
Fig. I5.5|) for the original log-normal coefficient configuration. The figure illustrates 
that successive approximations offer an accurate counterpart to the benchmark basis 
and perturbation values. 

5.2. Deterministic elliptic solution. In this subsection we assess the con- 
vergence of the respective solutions of Eq. (|2.1[) . In particular, we are interested in 
comparisons between the standard MsFEM solution u;, € Vh, and the proposed Ms- 
FEM solution uj^ € Vj,h- We first recall that Thm. 13.21 suggests convergence of 
to Uh as J — > oo. In order to verify this theoretical result we test two separate 
permeability configurations. In particular, we use a KL expansion (with a set of fixed 
random parameters) , and a single realization of a log- normal field analogous to those 
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Fig. 5.9. Fine scale KLE (left) and log-normal (right) coefficients posed on a 120 X 120 mesh 
(log scale) 



in Subsection 15. II However, the fine scale fields are now posed on at least a 120 x 120 
fine mesh. See Fig. 15.91 for a log-plot of each respective permeabililty field. Until oth- 
erwise noted, all MsFEM solutions in this subsection are obtained by using a 12 x 12 
coarse mesh. 

First, we consider a KLE coefficient which is posed on 120 x 120, 240 x 240, and 
360 x 360 fine meshes. For the variance we use a 1 = 2.25 and for the correlation 
lengths we use l x = 0.7 and l y = 0.04, respectively (see Eq. (|5.ip). For all examples 
we truncate the original KL expansion at n = 20 terms to obtain the full coefficient 
k = ko + ki . In order to test the convergence properties in Thm. I3.2l we are interested 
in computing a variety of proposed MsFEM solutions uj^ (for J = 1, . . .). Then, we 
compute the energy norm \\\uh — uj,h\\ | and the associated bound B u := (c) 1 / 2 (rj J+l + 

i1 r 7+r ) 1 2 |IM1I- We note that for |||u|||, we compute the energy norm of the standard 
fe'm elliptic solution. Fig. 15.101 illustrates the energy norm values and bounds for 
increasing J values, and two different KLE splitting configurations (m = 16 and 
m = 18). We note that, as before, a smaller r\ value yields a tighter bound. In addition, 
we see the pronounced convergence of the proposed MsFEM solution. Leaving the 
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Fig. 5.10. Energy error and error bound computations for the fine scale KLE coefficient; 
m = 16 terms (left), m = 18 terms (right); 120 X 120, 240 X 240, and 360 X 360 fine meshes 



Reference J = 0, Error = 3.6% J = 1, Error = 1.5% J = 2, Error = 0.8% 




Fig. 5.11. Convergence illustration for the reference MsFEM elliptic solution u^ (labeled 'Ref- 
erence'), and the corresponding new MsFEM elliptic solutions ujj x for J = 0, 1, 2; KLE coefficient, 
m = 12 



coarse mesh intact, Fig. 15.101 also serves to illustrate the effect that different local 
meshes have on the construction of the operators Mq and Mi from Thm. 12.11 In 
essence, the fine meshes give three successive local mesh refinements (from 10 x 10 to 
20 x 20 to 30 x 30) on which the local operators will be constructed. As shown in the 
figure, there is a negligible difference in the respective errors, thus illustrating that 
the fine scale on which the operators are computed does not have a significant effect 
on the error resulting from the proposed method. Furthermore, the results show that 
the original 120 x 120 mesh is sufficiently fine to resolve the necessary scales in the KL 
expansion. To further illustrate the convergence of the proposed MsFEM approach, 
we offer various elliptic solution plots in Fig. 15.111 In addition to the solution plots we 
offer the relative errors \ \\uh — uj^h\\\/\\\uh\\\- As J increases, we see that the relative 
error steadily decreases from 3.6% to 0.8%. 

We also consider the log-normal coefficient to test the convergence properties of 
the proposed approach. The energy norm and bound plots may be seen in Fig. 15.121 
where strength values of s c = 0.88 and s c = 0.9 are used. As expected, we see that the 
solutions converge and that a smaller 77 value yields a tighter bound. For this case, 
we also offer an illustration of various solution plots along with the relative errors 
values. See Fig. I5.13l for the convergence illustration. We note that as J, the relative 
error decreases from 1.0% to 0.2%. In this case we see that the initial error is less 
than its counterpart from the KLE coefficient. This is not unexpected since the KLE 
expansion represents a stronger form of heterogeneity (refer back to Fig. 15. 9j) . 

In addition to the results presented above, we also consider a variety of coarse 
mesh configurations for comparison. To recall, all previous examples in this subscc- 
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Fig. 5.12. Energy norm and error bound computations for the fine scale log-normal coefficient; 
■- 0.88 (left), s c = 0.9 (right) 
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Fig. 5.13. Convergence illustration for the reference MsFEM elliptic solution u^ (labeled 'Ref- 
erence'), and the corresponding new MsFEM elliptic solutions uj ^ for J = 0,1,2; log-normal 
coefficient, s c = 0.86 



tion were obtained from a fixed 12 x 12 coarse mesh. However, it is also fitting to 
illustrate the effects of different mesh configurations on the solution error. In particu- 
lar, Fig. 15.141 contains the errors quantities — Mj,/t||| (left) and \\\uh — uj,h\\\ (right) 
for a variety of mesh configurations and J values. The errors are obtained from the 
same KLE coefficient data with m = 16. We note that the left set of errors (the errors 
between the standard FEM solution and the proposed MsFEM solution) are essen- 
tially constant regardless of the value of J. In other words, the dominant source of 
error is clearly from the multiscale solution method, and the error between standard 
MsFEM and the proposed method is negligible. In addition, we see that a refinement 
of the coarse mesh yields a smaller error, which is to be expected from a multiscale 
solution technique. Of particular interest are the right set of errors \\\uh — uj,h\\\> 
which are computed using the same coarse mesh configurations. Most importantly, 
we see that a successive refinement of the coarse mesh does not significantly affect 
the error between standard MsFEM and the proposed method. We note that a re- 
fined coarse mesh does lead to a slight decrease in the error, however, it is clear from 
Fig. 15.141 that the proposed method is not sensitive with respect to the coarse mesh 
configuration. In particular, it is evident that J is the parameter which dictates the 
respective errors. 

5.3. Stochastic elliptic solution using the Monte Carlo method. In this 
subsection we address the stochastic problem described in Sect. |4] To begin, we are 
first interested in testing the stochastic bounds which are proved in Thm. 14.11 We 
remark that this stochastic result is analogous to the deterministic bounds offered in 



A splitting technique with applications to MsFEMs 



23 




Fig. 5.14. Energy error computations for the fine scale KLE coefficient; m = 16; 4 X 4, 12 X 
12, 20 X 20, and 30 X 30 coarse meshes 



Thm. 13.21 In order to compute the left hand side of the inequality in Thm. 14.11 we 
use the Monte Carlo approximation 

N N 

\\VkV(u h -uj, h )\\ L 2 (D) ® L 2 {a) « Y, W^«- u jm)\\lhd)/N = Y \\\<-Uj, h \\\/N, 

8=1 8=1 

(5.4) 

where the index s denotes a fixed sample value. In particular, we generate a stochastic 
field, compute the corresponding energy norms for s = 1, . . . , N, and obtain an average 
over the samples. For the right hand side of the inequality B u := (c) 1//2 (r? J+1 + 

ij~~7+i ) 1 ^ 2 \\\ / ^k^ u \\L 2 (D)<g>L 2 (n)i wc use the same type of Monte Carlo approximation 
for the stochastic integrals and we use the fully resolved FEM solution for u. In 
order to verify the bounds, we focus on a stochastic field which is generated from 
a truncated KL expansion. For the results in this subsection we employ correlation 
lengths of l x = l y = 0.1 and a variance of a 2 = 1.0 for the field construction. The 
series expansion is truncated at n = 20 terms, and the coefficient is posed on a 64 x 64 
fine mesh. We assume that the random coefficients (i = 1, . . . , 20) represent a 20 
dimensional vector in the hypercube [—1, l] 20 . In other words, for sampling we draw 
20 i.i.d. uniform random variables from the interval [—1, 1]. See Fig. l5.15l for a typical 
coefficient sample in this context. For all stochastic computations N — 200 samples 
are used, and for all MsFEM computations we use a 16 x 16 coarse mesh. 

To verify the bounds in Thm. 14.11 we offer two representative plots in Fig. 15.161 
The left hand side was obtained by keeping to = 14 terms in the original 20 term 
expansion, and the right hand side was obtained by keeping to = 18 terms in the orig- 
inal expansion. We note that more terms in the KL expansion (resulting in a smaller 
average -q value) yields a tighter bound, and that the expected energy norms deplete 
rapidly. These results are consistent with those obtained from the deterministic fields. 

To further illustrate the behavior of the stochastic problem we recall the energy 
ratio E(m) = ~= 1 ^)i-' , where Xi are the eigenvalues from the KL expansion. Wc 

note that as to increases this value is expected to quickly increase (in other words, 
1 — E(m) will quickly decrease). As more terms in the KL expansion typically yield 
smaller errors from the series approximation, we expect consistent behavior between 
1 — E{m) vs. the mean and variance quantities of the energy norm of the error, 
1 1 1 Uh — u j t h 1 1 1 • Fig. 15.171 illustrates the relationship between the statistics of the energy 
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Fig. 5.15. Typical elliptic coefficient sample from the KL expansion posed on a 64 X 64 mesh 
(log scale) 



m = 14 



m = 18 




Fig. 5.16. Expected energy norm and error bound computations for the KLE coefficient; m = 14 
terms (left), m = 18 terms (right) 



norm quantities and the energy ratio. For both plots we use J = 2 for the basis 
function approximations. For both the energy norm mean and variance comparisons, 
we note that as 1 — E(m) increases, the respective statistical values also increase. 
In particular, we see that less terms in the KL expansion yield errors that grow 
algebraically with respect to a decreasing energy ratio. 

To finish this subsection we offer statistical comparisons obtained from the ref- 
erence (standard) MsFEM solution Uh and the proposed MsFEM solution uj t h- See 
Fig. 15.181 for a comparison between the mean and variance of the respective elliptic 
solutions. We note that any differences are nearly undetectable, further verifying the 
accuracy of the proposed method. 

5.4. Stochastic elliptic solution using the parameter reduction colloca- 
tion. In this subsection we address the alternative parameter reduction collocation 
approach as described in Subsection 14.21 In particular, we are interested in testing 
the accuracy of Monte Carlo sampling with standard MsFEM, versus random pa- 
rameter reduction collocation with the proposed MsFEM approach using the Green's 
kernel. Throughout this section we recall that the total error can be decomposed into 
two main components. Namely, the error can be decomposed into the splitting error 
e sp i, and the collocation error e co ; = e co /(J, m,L). In order to assess the significance 
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Fig. 5.17. 1 — E(m) vs. the mean of ||| • | (left), and 1 — E(m) vs. the variance of 
(right); J = 2, Monte Carlo 
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Fig. 5.18. Statistical comparisons between standard MsFEM and the proposed MsFEM approach 



of the error contributions, we thoroughly test a variety of scenarios resulting from 
values of J (number of terms in the series expansion), to (splitting configuration), 
and L (collocation level). In a parallel setting the computational cost is decreased 
when the parameter reduction approach is implemented, and demonstrated accuracy 
of the technique would solidify it as a suitable sampling alternative. Throughout this 
subsection we use the same KLE configurations as in Subsection 15.31 

To motivate further disussion, we first offer Table 15.11 which compares the Monte 
Carlo sampling approach with the new parameter reduction sampling approach. The 
values in the table result from keeping to = 16 terms for construction of ko combined 
with interpolation levels L = 1, L = 2, and L = 3. In particular, we tabulate 
the values ^rc^n'r^ x 100% for 5 fixed 6 (random parameter) values, where Uh 
denotes the standard MsFEM solution with Monte Carlo sampling, and uj t h denotes 
the proposed MsFEM solution with the parameter reduction approach. As we can 
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see from Table 15.11 the parameter reduction approach closely recovers each individual 
sample of the elliptic solution, and a higher interpolation level leads to smaller errors. 
In all cases we note that the relative errors do not exceed 1%. 



Table 5.1 

Comparison of elliptic solution difference 



Q Sample 


Relative Errors (%) 


Level 1 


Level 2 


Level 3 


1 


0.64 


0.12 


0.05 


2 


0.30 


0.11 


0.10 


3 


0.33 


0.18 


0.17 


4 


0.67 


0.19 


0.08 


5 


0.25 


0.13 


0.12 



Although the individual sample results are promising, we are most interested in 
the statistical behavior of the respective solutions. In Fig. 15.191 we offer analogous 
results to those found in the previous subsection. Namely, we plot 1 — E(m) vs. 
the mean and variance quantities of the relative errors, \\\uh — ^J,h\\\/\\\uh\\\ and 
1 1 1 xt/i — uj^III/IIIm^IH, where we recall that Uj^ denotes the standard MsFEM solution 
technique combined with Monte Carlo sampling. In essence, we use these expected 
errors as benchmarks for comparison with the proposed method. More specifically, 
we are interested in comparing standard MsFEM with Monte Carlo sampling against 
the proposed MsFEM approach with parameter reduction collocation. We first note 
that an increase in interpolation level clearly yields a decrease in the relative error 
values. The level 2 interpolation errors nearly match the Monte Carlo results, and this 
slight discrepancy may be viewed as a trade of the increased efficiency of the proposed 
method. In other words, for a moderate collocation level it is natural to expect some 
minimal error contributions from collocation. However, the level 3 interpolation yields 
errors that are nearly identical to the Monte Carlo results. This is due to the fact 
that the collocation error is essentially dimished, and the only remaining error is due 
to the splitting. This will be discussed in further detail below. In Fig. 15.201 we single 
out the L = 3 results and plot them with the Monte Carlo results. We again note that 
the discrepancies are neglible, and point out the similarities in the statistical behavior 
which is illustrated in Fig. 15.171 In particular, we again encounter an increase of 
the mean and variance quantities with respect 1 — E(m) which is solely due to the 
splitting configuration. 

In addition to the high level of accuracy we obtain for a larger interpolation level, 
we emphasize that all plots from Fig. 15.191 illustrate relative errors that do not exceed 
1%. Even though the level 1 results do not closely match the Monte Carlo results, 
a neglible error of < 1% may be completely acceptable for many applications. In 
particular, the small vertical scale should be duely noted. From the L — 1 results we 
conclude that the proposed sampling method is not sensitive with respect to E(m), 
and fewer terms may be kept in the KL expansion without a significant loss in ac- 
curacy. This is due to the fact that e co i is dominant for the basic collocation level. 
These results can be viewed as a potential limitation and advantage of the approach 
if a level 1 interpolant is used. It may be a potential limitation because the addition 
of terms does not decrease the expected error significantly. However, these results 
may also be viewed as an advantage since 14 terms in a 20 term KL expansion (for 
example) yields similar errors as 18 terms in a 20 terms KL expansion. Thus, fewer 
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Fig. 5.19. 1 — E(m) vs. the expected value (left), and variance (right) of the relative error 
quantities; J = I, L = 1,2,3 and Monte Carlo 
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Fig. 5.20. 1 — E(m) vs. the expected value (left), and variance (right) of the relative error 
quantities; J = 1, L = 3 and Monte Carlo 



terms may be used in the decomposition with a ncglible loss in accuracy. However, as 
the results show, no sacrifice in accuracy is necessary if a higher interpolation level is 
used. This would simply amount to more precomputation steps in which additional 
sparse grid points are considered for the Smolyak intcrpolant. 

We also test the sensitivity of the approach with respect to the number of terms 
J which are kept in Eq. ()4.6|) . All previous results in this subsection were obtained 
from a value of J = 1 in the expansions; however, it is fitting to offer a comparison 
between solutions obtained by keeping more terms in the expansion. For the following 
comparisons we use the same level 1 interpolation results as seen in Fig. 15.191 and test 
the results corresponding to J = 1 and J = 4. In Fig. 15.211 we note that the method 
does not exhibit sensitivity with respect to J. This is again due to the fact that 
the collocation error is dominant for a level L = 1 interpolation. We see that the 
mean errors may be slightly decreased by keeping more terms in the series expansion, 
yet this increase in accuracy is subtle. This may be attributed to the fact that 
the L = 1 Green's function intcrpolant I m G(&o) in Eq. (|4.6[) does not guarantee 
faster convergence depending on the number of terms kept in the series expansion. 
More specifically, the dominant collocation error is inhereted through the iterative 
procedure. As more terms do not yield a significant gain in accuracy, keeping less 
terms is preferable in this low interpolation level setting. However, as the splitting 
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Fig. 5.21. Comparison between the expected value (left), and variance (right) of the relative 
error quantities; L = 1, J = 1, 4 



error becomes dominant (i.e., a higher interpolation level is considered) using more J 
terms would, in fact, result in more pronounced accuracy. We next elaborate on some 
additional effects of a higher interpolation level. 

To finish this section we offer detailed relative error comparisons for \\\uh — 
uj,h\\\/\\\uh\\ \ (total error), (splitting error) , and 

(collocation error). In doing so, the aim is to solidify the contention that a higher 
interpolation level indeed diminishes the the collocation error. In turn, if the neglible 
total error (< 1%) which results from a lower level interpolation (e.g., L = 1) is not 
suitable, a higher level interpolation (e.g., L = 3) may be used to negate the colloca- 
tion error. To reiterate, uj.h denotes a solution obtained from the proposed MsFEM 
approach with Monte Carlo sampling, and uj.h denotes a solution obtained from the 
proposed MsFEM approach with parameter reduction. See Fig. 15.221 for an illustra- 
tion of the various errors. We introduce a slight abuse of the original notation from 
Subsection 14.21 and use e to denote the total relative error, e sp i to denote the relative 
splitting error, and e co ; to denote the relative collocation error in Fig. 15.221 In the 
left hand side of Fig. 15.221 we note that the smallest errors result from the proposed 
MsFEM combined with Monte Carlo sampling (i.e., the splitting error is small). In 
addition, the total error and collocation error are comparable. Thus, we conclude that 
the collocation approach yields the dominant soure of error for a low level interpolant. 
However, we note the significant difference in the right hand side of Fig. 15.221 In this 
case we note two important factors. First, the total error from a L = 3 interpolant 
is smaller than its L = 1 counterpart (as expected). Furthermore, we see that the 
splitting error is now the dominant source of error, and the collocation error is much 
smaller. In other words, an increase in the interpolation level accomplishes the task 
of significantly reducing the collocation error. Thus, we conclude that for a higher 
interpolation level the parameter reduction approach behaves much like the standard 
Monte Carlo counterpart due to the minimal effect of collocation error. 

6. Conclusions. In this paper we present a new MsFEM approach for solving 
elliptic equations with coefficients that vary on many length scales and contain uncer- 
tainties. Through considering a coefficient decomposition combined with a Green's 
function approach, we are able to construct new MsFEM solutions which closely re- 
cover traditional MsFEM solutions. In a deterministic setting, rigorous error estimates 
and bounds are first presented for a representative basis function within the MsFEM 
approximation space. Using the initial basis function results, we offer a rigorous error 



A splitting technique with applications to MsFEMs 29 



Error Contributions [L = 1) Error Contributions [L — 3) 




Sample Sample 

Fig. 5.22. Comparison between the relative error quantities resulting from respective methods; 
m = 16, J = 1, L = 1 (left), L = 3 (right) 

analysis describing the behavior of the elliptic solutions that are sought in the space 
that is spanned by the new multiscalc basis functions. Under appropriate assumptions 
on the coefficient splitting configuration, we are ultimately able to construct approxi- 
mate solutions which nicely converge to a benchmark solution. The basis function and 
elliptic solution analysis are thoroughly verified through a number of representative 
numerical examples. In a stochastic setting, the proposed solution method is shown 
to reduce the number of sample solutions that must be constructed for assessing the 
statistical behavior of the system. In particular, the splitting gives rise to a situa- 
tion where the random parameter dimension can be reduced, and where stochastic 
collocation becomes an efficient alternative to direct Monte Carlo sampling. As the 
parameter reduction sampling approach involves a number of pre-computation steps 
that are completely independent, this approach is especially desirable in a parallel 
setting. Analogous error bounds are derived for the stochastic problem, and the suc- 
cessful performance of the proposed method is verified through a variety of numerical 
examples. 
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Appendix A. proof of Lemma 13.31 Let £o solve equation (|2.6[) . Then inte- 
gration by parts gives 

(*oV|o,v|o) = (fciV(n-/)z,v£ ) = (^ v(n-/)z,vfo). 

This implies that 

||VfeuVeo||o,K <^|lv / ^V(n-/);|| , /f <2r, A -||v / ^V/||o,K. (A.l) 
Let £j (j = 1, 2, ■ • • ) solve Eq. (|2.7j) . Then a similar argument shows 

\\Vk^Vi 3 \\o,K < VK\\VkoVij-i\\o,K- (A.2) 
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Using (|A.2|) recursively and (|A.1[) . we have 

\\VkoVii\\o,K < Tk-IIV^V&Ho.jc < 2?7 A - +1 || VW|| .a, j = 0,1,2,- •• . 
The proof is completed. 

Appendix B. proof of Lemma 13.41 Since <fij = (I — H)l + £,/ and 

(7 — n)Z + £, it suffices to show that 

"II* = 0. 

and the sequence of equations in (|2.7p , it follows that 



lim IIIO 

J— f oo 



By adding Eq. 



j-i 



-V • 



= v • fciV(X)li) - v • (fciV(n- 1)0 in # 



(B.l) 



£6 =0 on ax. 



Because J2j=o £? = O — £j> ^l- (EH is simplified to 

J -V • (fcV£j) = -V • (fciVO) - V • (fciV(n - 1)1) in K 
\ O = on dK. 

Applying integration by parts to (|B.2|) . we have 

(fcvo, vo) - (ftiVfj, vo) + (fciV(n - /)/, vf j) 



(B.2) 



fe v(n - j)/, Vfeve 



As a consequence, it follows immediately that 
IINIk< 



IU-(/f)llv / ^ve/||o,7< + ||^=||L-(A-)llv / ^v(n-/);||o,K 



<2|| 



'fcfco 

fci 
s/kko 
h 



'kko 



IU~(K)(< +1 + 1)||VW|| , 



.A' 



IU~(a-)||v fcoVZ||o,jr as J -> oo, 



where we have used Lemma 13.31 in the last step. Hence, £j is convergent, and <f>j is 
convergent as well. Subtracting Eq. (|B.2|) from Eq. (|2.5[) . we have 



-V • (fcV(£ - O)) = V • (fcVf j) in /v 
£ - £j = on dK. 



(B.3) 



Performing integration by parts and using the Cauchy-Schwarz inequality for Eq. 
(IB. 31) . then we obtain 



me — e 



J K 



< 



(A-)||\AoV£/||o. 



K 



<2||-fe=|| i - (J0 ^+ 1 | 



as J -> 0. 



%V/|| 



0,K 
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This completes the proof. 

Appendix C. proof of Proposition 13.51 Because <fij = (I — H)l + £j and 

(j> = (7 — n)^ + £, it suffices to show 



\\\t-tj\\\K<c l ^ n j+*. 



In fact, the proof of Lemma 13^1 implies that 



Ollk<2|| 



ki 



(K)ri J K +1 \\\/koVl\\o,K 



o 



<2|||^|Uoo W ^||V^| 



l- (jo ||V/||o,jc = Ci^kv J K +2 



This completes the proof. 
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